#Ejercicios {.tabset .tabset-fade .tabset-pills}
###Punto 10
(a) Produzaca algunos resumenes numéricos y gráficos de los datos Weekly.¿ parace que hay patrones?.
Resumen de variables:
library(ISLR)
data("Weekly")
summary(Weekly)
## Year Lag1 Lag2 Lag3
## Min. :1990 Min. :-18.1950 Min. :-18.1950 Min. :-18.1950
## 1st Qu.:1995 1st Qu.: -1.1540 1st Qu.: -1.1540 1st Qu.: -1.1580
## Median :2000 Median : 0.2410 Median : 0.2410 Median : 0.2410
## Mean :2000 Mean : 0.1506 Mean : 0.1511 Mean : 0.1472
## 3rd Qu.:2005 3rd Qu.: 1.4050 3rd Qu.: 1.4090 3rd Qu.: 1.4090
## Max. :2010 Max. : 12.0260 Max. : 12.0260 Max. : 12.0260
## Lag4 Lag5 Volume
## Min. :-18.1950 Min. :-18.1950 Min. :0.08747
## 1st Qu.: -1.1580 1st Qu.: -1.1660 1st Qu.:0.33202
## Median : 0.2380 Median : 0.2340 Median :1.00268
## Mean : 0.1458 Mean : 0.1399 Mean :1.57462
## 3rd Qu.: 1.4090 3rd Qu.: 1.4050 3rd Qu.:2.05373
## Max. : 12.0260 Max. : 12.0260 Max. :9.32821
## Today Direction
## Min. :-18.1950 Down:484
## 1st Qu.: -1.1540 Up :605
## Median : 0.2410
## Mean : 0.1499
## 3rd Qu.: 1.4050
## Max. : 12.0260
Matriz de gráfico de dispersión:
pairs(Weekly)
En la primera matriz muestra un resumen de las variables a trabajar, observamos que todas estan a una misma escala, ecepto la varible dirección que cuenta con dos niveles (Alto y bajo). En la segunda matriz las varibles lag1,lag2,lag3,lag4,lag5 y volumen, son variables que entre ellas aparentemente las observaciones se concentrar en el centro; podemos notar que las varible año y variable volume tiene un comportamiento de un polinomio de grado dos acsendente mientras pasa el tiempo.
Matriz de correlaciones
cor(Weekly[,-9])
## Year Lag1 Lag2 Lag3 Lag4
## Year 1.00000000 -0.032289274 -0.03339001 -0.03000649 -0.031127923
## Lag1 -0.03228927 1.000000000 -0.07485305 0.05863568 -0.071273876
## Lag2 -0.03339001 -0.074853051 1.00000000 -0.07572091 0.058381535
## Lag3 -0.03000649 0.058635682 -0.07572091 1.00000000 -0.075395865
## Lag4 -0.03112792 -0.071273876 0.05838153 -0.07539587 1.000000000
## Lag5 -0.03051910 -0.008183096 -0.07249948 0.06065717 -0.075675027
## Volume 0.84194162 -0.064951313 -0.08551314 -0.06928771 -0.061074617
## Today -0.03245989 -0.075031842 0.05916672 -0.07124364 -0.007825873
## Lag5 Volume Today
## Year -0.030519101 0.84194162 -0.032459894
## Lag1 -0.008183096 -0.06495131 -0.075031842
## Lag2 -0.072499482 -0.08551314 0.059166717
## Lag3 0.060657175 -0.06928771 -0.071243639
## Lag4 -0.075675027 -0.06107462 -0.007825873
## Lag5 1.000000000 -0.05851741 0.011012698
## Volume -0.058517414 1.00000000 -0.033077783
## Today 0.011012698 -0.03307778 1.000000000
(b) Utilice el conjunto de datos completo para realizar una regresion logística con Direction como la respuesta y la cinco varibles de retraso mas Volume como predictoras. Use la funcion summary para imprimir los resultados. ¿algunos de los predictores parece ser estadisticamente significativos? si es asi ¿cuales son estas variables?
Resumen Regresión logística
mod1 <- glm( Direction ~ Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data = Weekly,family = "binomial")
summary(mod1)
##
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 +
## Volume, family = "binomial", data = Weekly)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.6949 -1.2565 0.9913 1.0849 1.4579
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.26686 0.08593 3.106 0.0019 **
## Lag1 -0.04127 0.02641 -1.563 0.1181
## Lag2 0.05844 0.02686 2.175 0.0296 *
## Lag3 -0.01606 0.02666 -0.602 0.5469
## Lag4 -0.02779 0.02646 -1.050 0.2937
## Lag5 -0.01447 0.02638 -0.549 0.5833
## Volume -0.02274 0.03690 -0.616 0.5377
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1496.2 on 1088 degrees of freedom
## Residual deviance: 1486.4 on 1082 degrees of freedom
## AIC: 1500.4
##
## Number of Fisher Scoring iterations: 4
Podemos observar que la variable Lag2 con un P-valor de 0.0296 es estadisiticamente significativa con un nivel de confianza de 0.05.
(c) Compute la matrix de confusión y la fracción general de las predicciones correctas. Explica qué te dice la matriz de confusión logística.
Tabla de confusión
prediction <- mod1$fitted.values
pred<- rep("Dow",length(prediction))
pred[prediction > 0.5]<- "Up"
table(pred,Weekly$Direction)
##
## pred Down Up
## Dow 54 48
## Up 430 557
Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{54+557}{1089} = 0.5610 \times 100\) = 56.10% del tiempo . La tasa de error con los datos de entrenamiento es del 100%-56.10% = 43.8%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{557}{48+557}=0.92\times100%\)=92% de la veces, tambien cuando el mercado baja el modelo acierta un \(\frac{54}{54+430}=0.1115\times 100\) = 11.15% del tiempo.
(d) Ahora ajuste el modelo de regresión logística utilizando un período de datos de entrenamiento de 1990 a 2008, con Lag2 como el unico predictor. Calcule la matriz de confusión y la fracción general de las predicciones correctas para los datos retenidos (Es decir, los datos de 2009 y 2010)
Resumen Modelo con lag2.
train <- subset.data.frame(x = Weekly,subset = Year < 2009)
test2009_2010 <- subset.data.frame(Weekly,subset = Year >=2009)
mod2 <- glm(Direction ~ Lag2 ,data = train,family = binomial)
summary(mod2)
##
## Call:
## glm(formula = Direction ~ Lag2, family = binomial, data = train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.536 -1.264 1.021 1.091 1.368
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.20326 0.06428 3.162 0.00157 **
## Lag2 0.05810 0.02870 2.024 0.04298 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 1354.7 on 984 degrees of freedom
## Residual deviance: 1350.5 on 983 degrees of freedom
## AIC: 1354.5
##
## Number of Fisher Scoring iterations: 4
Tabla de confusión
prediction2 <- predict(object = mod2,newdata = test2009_2010,type = "response")
pred2<- rep("Dow",length(prediction2))
pred2[prediction2 > 0.5]<- "Up"
table(pred2,test2009_2010$Direction)
##
## pred2 Down Up
## Dow 9 5
## Up 34 56
Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{9+56}{104}=0.625\times 100\)= 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-62.5% = 37.5%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{56}{5+56}=0.9180\times100%\)=91.8% del tiempo, tambien cuando el mercado baja el modelo acierta un \(\frac{34}{34+9}=0.79\times 100\) = 79.6% del tiempo.
(e) Repita (d) usando LDA
Resumen modelo LDA:
library(MASS)
modlda <- lda(Direction~Lag2,data = train)
modlda
## Call:
## lda(Direction ~ Lag2, data = train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
##
## Coefficients of linear discriminants:
## LD1
## Lag2 0.4414162
Tabla de confusión:
prediction_LDA <- predict(object = modlda,newdata = test2009_2010)
table(prediction_LDA$class,test2009_2010$Direction)
##
## Down Up
## Down 9 5
## Up 34 56
Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{9+56}{104}=0.625\times 100\)= 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-62.5% = 37.5%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{56}{5+56}=0.9180\times100%\)=91.8% de la veces, tambien cuando el mercado baja el modelo acierta un \(\frac{34}{34+9}=0.79\times 100\) = 79.6% del tiempo.
(f) Repita (d) usando QDA
Resumen modelo QDA
mod_QDA <- qda(Direction~Lag2,data = train)
mod_QDA
## Call:
## qda(Direction ~ Lag2, data = train)
##
## Prior probabilities of groups:
## Down Up
## 0.4477157 0.5522843
##
## Group means:
## Lag2
## Down -0.03568254
## Up 0.26036581
prediction_QDA <- predict(mod_QDA,newdata = test2009_2010)
table(prediction_QDA$class,test2009_2010$Direction)
##
## Down Up
## Down 0 0
## Up 43 61
Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de \(\frac{0+61}{104}=0.5865 \times 100\)= 58.65% del tiempo. La tasa de error con los datos de prueba es del 100%-58.65% = 41.35%. Tambien apreciamos que cuando el mercado sube, El modelo acierta \(\frac{61}{61}=1\) osea el 100% de la veces, Se resalta este modelo de anialisis de discriminación cudratica elige una dirección alta todo el tiempo.
(g) Repita (d) usando Knn con k=1.
Table de confusión:
library(class)
train.x <- as.matrix(train$Lag2)
test.x <- as.matrix(test2009_2010$Lag2)
train.Direction <- train$Direction
set.seed(1)
prediction_knn <- knn(train = train.x,test = test.x,cl = train.Direction,k = 1)
table(prediction_knn,test2009_2010$Direction)
##
## prediction_knn Down Up
## Down 21 30
## Up 22 31
Con esta tabla de confusión y los datos de pruba, podemos concluir que el porcentaje de acierto en la prediccion es de 62.5% del tiempo. La tasa de error con los datos de prueba es del 100%-56.10% = 43.9%. Tambien apreciamos que cuando el mercado sube, El modelo acierta 50.81% de la veces, tambien cuando el mercado baja el modelo acierta un 48.83% del tiempo.
(h) ¿ cuál de estos métodos parece proporcionar los mejores resultados con estos datos ?
Con los modelos anteriores, vemos que el modelo de regresión logístico y LDA, La tasa de error es mínima. tambien para los modelos QDA y KNN un poco menores.
(i) Experimente con diferentes combinaciones de predictores, incluyendo posibles transformaciones, para cada uno de los métodos, reporte la variables,método y la matriz de confusión asociada que parece proporcionar los mejores resultados en los mejores resultados en los datos retenidos.Tenga encuenta también debe experimentar con los valores para K en la clasificación con KNN.
La variable mas significativa es Lag2 como se mostro en los anteriores puntos para tener una interacción con las segunda varible significativa.
library(MASS)
library(class)
mod_inter_1 <- glm(Direction~Lag2:Lag1,family = binomial,data = train)
prediction_inter_1 <- predict(mod_inter_1,test2009_2010,type = "response")
pred_inter_1 <- rep("Dow",length(prediction_inter_1))
pred_inter_1[prediction_inter_1 > 0.5] <- "Up"
mod_lda_inter <-lda(Direction~Lag2:Lag1,data = train)
prediction_LDA_inter <- predict(mod_lda_inter,test2009_2010)
mod_QDA_inter <- qda(Direction ~ Lag2+sqrt(abs(Lag2)),train)
pred_QDA_inter <- predict(mod_QDA_inter,test2009_2010)
set.seed(0511)
#knn k=10
prediction_knn1 <- knn(train = train.x,test = test.x,cl = train.Direction,k = 10)
#knn=100
prediction_knn2 <- knn(train = train.x,test = test.x,cl = train.Direction,k = 100)
Tablas de modelos realizados con interacciones
| Modelo | predicción correcta | tasa de presición cuando el mercado aumenta |
|---|---|---|
| regresión logistica con interacción | 57.69% | 98.36% |
| LDA con interacción | 57.69% | 98.36%% |
| QDA con \(\sqrt(abs(Lag2))\) | 57.69% | 78.68% |
| knn con K=10 | 57.69% | 68.85% |
| Knn con k=100 | 56.73% | 80.32% |
De esta tabla se puede concluir que el modelo de regresión logística y LDA tienen el mejor rendimiento en tasas de error de prueba.
(a) Cree una variable binaria, mpg01, que contenga un 1 si mpg contiene un valor por encima de su mediana, y un 0 si mpg contiene un valor por debajo de su mediana. Puede calcular la mediana utilizando la función median(). Tenga en cuenta que puede resultarle útil utilizar la función data.frame() para crear un único conjunto de datos que contenga tanto mpg01 como las otras variables automáticas.
library(ISLR)
data(Auto)
n <- length(Auto$mpg)
mpg01 <- ifelse( Auto$mpg > median(Auto$mpg),1,0)
Auto <- data.frame(Auto,mpg01)
Auto$mpg01 <- as.factor(Auto$mpg01)
knitr::kable(head(Auto),caption = "Cabeza de base de datos Auto con variable mpg01")
| mpg | cylinders | displacement | horsepower | weight | acceleration | year | origin | name | mpg01 |
|---|---|---|---|---|---|---|---|---|---|
| 18 | 8 | 307 | 130 | 3504 | 12.0 | 70 | 1 | chevrolet chevelle malibu | 0 |
| 15 | 8 | 350 | 165 | 3693 | 11.5 | 70 | 1 | buick skylark 320 | 0 |
| 18 | 8 | 318 | 150 | 3436 | 11.0 | 70 | 1 | plymouth satellite | 0 |
| 16 | 8 | 304 | 150 | 3433 | 12.0 | 70 | 1 | amc rebel sst | 0 |
| 17 | 8 | 302 | 140 | 3449 | 10.5 | 70 | 1 | ford torino | 0 |
| 15 | 8 | 429 | 198 | 4341 | 10.0 | 70 | 1 | ford galaxie 500 | 0 |
(b) Explore los datos gráficamente para investigar la asociación entre mpg01 y las otras características. ¿Cuál de las otras características parece más útil para predecir mpg01? Los diagramas de dispersión y los diagramas de caja pueden ser herramientas útiles para responder a esta pregunta. Describe tus hallazgos.
Correlaciones entre variables:
Auto1 <- Auto
Auto1$mpg01 <- as.numeric(Auto1$mpg01)
cor(Auto1[,-9])
## mpg cylinders displacement horsepower weight
## mpg 1.0000000 -0.7776175 -0.8051269 -0.7784268 -0.8322442
## cylinders -0.7776175 1.0000000 0.9508233 0.8429834 0.8975273
## displacement -0.8051269 0.9508233 1.0000000 0.8972570 0.9329944
## horsepower -0.7784268 0.8429834 0.8972570 1.0000000 0.8645377
## weight -0.8322442 0.8975273 0.9329944 0.8645377 1.0000000
## acceleration 0.4233285 -0.5046834 -0.5438005 -0.6891955 -0.4168392
## year 0.5805410 -0.3456474 -0.3698552 -0.4163615 -0.3091199
## origin 0.5652088 -0.5689316 -0.6145351 -0.4551715 -0.5850054
## mpg01 0.8369392 -0.7591939 -0.7534766 -0.6670526 -0.7577566
## acceleration year origin mpg01
## mpg 0.4233285 0.5805410 0.5652088 0.8369392
## cylinders -0.5046834 -0.3456474 -0.5689316 -0.7591939
## displacement -0.5438005 -0.3698552 -0.6145351 -0.7534766
## horsepower -0.6891955 -0.4163615 -0.4551715 -0.6670526
## weight -0.4168392 -0.3091199 -0.5850054 -0.7577566
## acceleration 1.0000000 0.2903161 0.2127458 0.3468215
## year 0.2903161 1.0000000 0.1815277 0.4299042
## origin 0.2127458 0.1815277 1.0000000 0.5136984
## mpg01 0.3468215 0.4299042 0.5136984 1.0000000
Grafico de dispersion
pairs(Auto)
attach(Auto)
## The following object is masked _by_ .GlobalEnv:
##
## mpg01
layout(rbind(c(1,1,1,2,2,2,3,3,3),c(4,4,4,5,5,5,6,6,6)))
boxplot(cylinders~mpg01,data = Auto,main="Cilindros vs mpg01")
boxplot(displacement ~mpg01,data = Auto,main="Desplazamiento vs mpg01")
boxplot(horsepower ~mpg01,data = Auto,main="Caballos de fuerza vs mpg01")
boxplot(weight ~mpg01,data = Auto,main="Peso del vehiculo vs mpg01")
boxplot(acceleration ~mpg01,data = Auto,main="Aceleración vs mpg01")
boxplot(year ~mpg01,data = Auto,main="Modelo vs mpg01")
(c) Divida los datos en un conjunto de entrenamiento y un conjunto de prueba.
train <- sample(seq(length(Auto$mpg01)),length(Auto$mpg01)*0.70,replace = FALSE)
data_train <- Auto[train,]
data_test <- Auto[-train,]
| Entranamiento | Prueba | Total |
|---|---|---|
| 274 | 118 | 392 |
| ————- | ——- | —– |
(d) Realice LDA en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?
require(MASS)
fit.lda <- lda(mpg01 ~ cylinders+displacement+horsepower+weight,Auto)
fit.lda
## Call:
## lda(mpg01 ~ cylinders + displacement + horsepower + weight, data = Auto)
##
## Prior probabilities of groups:
## 0 1
## 0.5 0.5
##
## Group means:
## cylinders displacement horsepower weight
## 0 6.765306 273.1582 130.11224 3620.403
## 1 4.178571 115.6658 78.82653 2334.765
##
## Coefficients of linear discriminants:
## LD1
## cylinders -0.4708126926
## displacement -0.0014339977
## horsepower 0.0044241441
## weight -0.0009846242
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$mpg01)
##
## 0 1
## 0 51 3
## 1 11 53
mean(pred.lda$class != data_test$mpg01)
## [1] 0.1186441
Con los resultados anteriores podemos concluir que la tasa de error es del 11.01%.
(e) Realice QDA en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?
fit.qda <- qda(mpg01~cylinders+displacement+horsepower,data = data_train)
fit.qda
## Call:
## qda(mpg01 ~ cylinders + displacement + horsepower, data = data_train)
##
## Prior probabilities of groups:
## 0 1
## 0.4890511 0.5109489
##
## Group means:
## cylinders displacement horsepower
## 0 6.783582 278.5970 132.8881
## 1 4.185714 117.8071 80.1500
pred.qda <- predict(fit.qda,data_test)
table(pred.qda$class,data_test$mpg01)
##
## 0 1
## 0 53 3
## 1 9 53
mean(pred.qda$class != data_test$mpg01)
## [1] 0.1016949
Con los resultados anteriores podemos concluir que la tasa de error es del 11.01%.
(f) Realice una regresión logística en los datos de entrenamiento para predecir mpg01 usando las variables que parecían más asociadas con mpg01 en (b). ¿Cuál es el error de prueba del modelo obtenido?
fit.glm <- glm(mpg01 ~ cylinders+displacement+horsepower+weight+acceleration,family = binomial,data_train)
summary(fit.glm)
##
## Call:
## glm(formula = mpg01 ~ cylinders + displacement + horsepower +
## weight + acceleration, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -2.4925 -0.1566 0.1188 0.3611 3.2309
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 11.730488 3.276787 3.580 0.000344 ***
## cylinders -0.136576 0.404783 -0.337 0.735811
## displacement -0.007845 0.009597 -0.817 0.413649
## horsepower -0.033194 0.024676 -1.345 0.178556
## weight -0.002468 0.001116 -2.212 0.026996 *
## acceleration 0.038322 0.159330 0.241 0.809930
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 379.71 on 273 degrees of freedom
## Residual deviance: 147.81 on 268 degrees of freedom
## AIC: 159.81
##
## Number of Fisher Scoring iterations: 7
pr<- predict(fit.glm,data_test,type = "response")
pred.glm <- rep(0,length(pr))
pred.glm[pr>0.5]<-1
table(pred.glm,data_test$mpg01)
##
## pred.glm 0 1
## 0 52 3
## 1 10 53
mean(pred.glm != data_test$mpg01)
## [1] 0.1101695
Con los resultados anteriores podemos concluir que la tasa de error es del 11.86%.
(g) Realice KNN en los datos de entrenamiento, con varios valores de K, para predecir mpg01. Use solo las variables que parecían más asociadas con mpg01 en (b). ¿Qué errores de prueba obtienes? ¿Qué valor de K parece tener el mejor rendimiento en este conjunto de datos?
require(ggplot2)
## Loading required package: ggplot2
##
## Attaching package: 'ggplot2'
## The following object is masked from 'Auto':
##
## mpg
require(factoextra)
## Loading required package: factoextra
## Welcome! Related Books: `Practical Guide To Cluster Analysis in R` at https://goo.gl/13EFCZ
require(class)
require(caret)
## Loading required package: caret
## Loading required package: lattice
set.seed(0511)
ctrl <- trainControl(method="repeatedcv",repeats = 3)
knnFit <- train(mpg01 ~ cylinders+displacement+horsepower+weight+acceleration, data =data_train[,-c(1,7,8,9)], method = "knn", trControl = ctrl, preProcess = c("center","scale"), tuneLength = 20)
knnFit
## k-Nearest Neighbors
##
## 274 samples
## 5 predictor
## 2 classes: '0', '1'
##
## Pre-processing: centered (5), scaled (5)
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 247, 247, 247, 247, 246, 246, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 5 0.9063492 0.8123034
## 7 0.9003086 0.8002292
## 9 0.8966931 0.7931661
## 11 0.9039683 0.8077058
## 13 0.9039683 0.8077058
## 15 0.9039683 0.8077058
## 17 0.9039683 0.8076792
## 19 0.9039683 0.8076792
## 21 0.9015432 0.8028595
## 23 0.9015432 0.8028595
## 25 0.9015432 0.8028595
## 27 0.9015432 0.8028595
## 29 0.9015432 0.8028595
## 31 0.9015432 0.8028595
## 33 0.9015432 0.8028595
## 35 0.9015432 0.8028595
## 37 0.9015432 0.8028595
## 39 0.9027337 0.8052404
## 41 0.9027337 0.8052404
## 43 0.9015432 0.8028595
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.
se obtuvo un k optimo de 7 procedemos con el modelo.
train.x <- as.matrix(data_train[,-c(1,7,8,9)])
test.x <- as.matrix(data_test[,-c(1,7,8,9)])
train.y <- as.matrix(data_train$mpg01)
test.y <- as.matrix(data_test$mpg01)
set.seed(0511)
pred.knn <- knn(train = train.x,test = test.x,cl = train.y,k = 7)
table(pred.knn,test.y)
## test.y
## pred.knn 0 1
## 0 52 4
## 1 10 52
mean(pred.knn !=test.y)
## [1] 0.1186441
con los resultados obtenidos con un K=7 optimo la tasa de error es de 13.55%.
Conclución general: Los mejores modelos son LDA y QDA con una tasa de error de 11.01%
(a) Escriba una función, Power(), que imprima el resultado de elevar 2 a la tercera potencia. En otras palabras, su función debe calcular \(2^3\) e imprimir los resultados.Sugerencia: recuerde que \(x^a\) eleva x a la potencia a. Use la función print() para generar el resultado.
Power <- function() {
2^3
}
Power()
## [1] 8
(b) Cree una nueva función, Power2(), que le permita pasar dos números, x y a, e imprima el valor de \(x^a\). Puede hacer esto comenzando su función con la línea
power <- function(){2^3}
power()
## [1] 8
(b) Cree una nueva función, Power2(), que le permita pasar dos números, x y a, e imprima el valor de \(x^a\). Puede hacer esto comenzando su función con la línea
Power2 <- function (x,a){x^a}
Power2(3,8)
## [1] 6561
Power2 <- function(x, a) {
x^a
}
Power2(3, 8)
## [1] 6561
(c) Usando la función Power2() que acaba de escribir, calcule \(10^3\), \(8^17\) y \(131^3\).
Power2(10,3)
## [1] 1000
Power2(8,17)
## [1] 2.2518e+15
Power2(131,3)
## [1] 2248091
(d) Ahora cree una nueva función, Power3(), que realmente devuelve el resultado x^a como un objeto R, en lugar de simplemente imprimirlo en la pantalla. Es decir, si almacena el valor x^a en un objeto llamado resultado dentro de su función, simplemente puede return() este resultado, utilizando la siguiente línea:
power3 <- function(x,a){
resultado <- x^a
return(resultado)
}
Power3 <- function(x , a) {
result <- x^a
return(result)
}
(e) Ahora, utilizando la función Power3 (), cree un plot de f(x) = \(x^2\). El eje x debe mostrar un rango de enteros del 1 al 10, y el eje y debería mostrar \(x^2\). Rotule los ejes apropiadamente y use un título apropiado para la figura. Considere mostrar el eje x, el eje y o ambos en la escala logarítmica. Puede hacerlo utilizando log = ‘‘ x’’, log = ‘‘ y’’ o log = ‘‘ xy’’ como argumentos de la función plot()
función a graficar:
\[f(x)=x^2\]
x <- 1:10
plot(x,power3(x,2),log = "xy",xlab = "log de x",ylab = "log de x^2",main = "log de x^2 vs log de x")
x <- 1:10
plot(x, Power3(x, 2), log = "xy", xlab = "Log of x", ylab = "Log of x^2", main = "Log of x^2 vs Log of x")
(f) Cree una función, PlotPower (), que le permite crear una gráfica de x contra x ^ a para un a fijo y para un rango de valores de x. Por ejemplo, si llamas
entonces se debe crear un gráfico con un eje x que tome valores 1,2, …, 10 y un eje y que toma los valores \(1^3\),\(2^3\), …, \(10^3\).
plotpower <- function(x,a){
plot(x,power3(x,a))
}
plotpower(1:10,3)
PlotPower <- function(x, a) {
plot(x, Power3(x, a))
}
PlotPower(1:10, 3)
Utilizando el conjunto de datos de Boston, ajuste los modelos de clasificación para predecir si un suburbio determinado tiene una tasa de criminalidad superior o inferior a la mediana. Explore los modelos de regresión logística, LDA y KNN utilizando varios subconjuntos de predictores. Describe tus hallazgos.
library(MASS)
data("Boston")
medianacrim <- median(Boston$crim)
crim01 <- rep(0,length(Boston$crim))
crim01[Boston[1]>medianacrim] <- 1
Boston$crim01 <- crim01
train <- sample(seq(length(Boston$crim01)),length(Boston$crim01)*0.70,replace = FALSE)
data_train <- Boston[train,]
data_test <- Boston[-train,]
Regresion logistica:
fit.glm <- glm(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat,family = binomial,data_train)
summary(fit.glm)
##
## Call:
## glm(formula = crim01 ~ zn + indus + chas + nox + rm + age + dis +
## rad + tax + ptratio + black + lstat, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.9848 -0.1797 -0.0017 0.0034 3.3317
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -30.570448 7.823691 -3.907 9.33e-05 ***
## zn -0.066360 0.039966 -1.660 0.096830 .
## indus -0.043417 0.053660 -0.809 0.418451
## chas 0.812084 0.893650 0.909 0.363495
## nox 46.839608 8.845996 5.295 1.19e-07 ***
## rm 0.280004 0.563847 0.497 0.619474
## age 0.022823 0.014516 1.572 0.115881
## dis 0.483227 0.236176 2.046 0.040752 *
## rad 0.693487 0.184343 3.762 0.000169 ***
## tax -0.009323 0.003294 -2.830 0.004649 **
## ptratio 0.178923 0.121506 1.473 0.140874
## black -0.006888 0.005393 -1.277 0.201544
## lstat -0.036400 0.057403 -0.634 0.526006
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 490.70 on 353 degrees of freedom
## Residual deviance: 151.08 on 341 degrees of freedom
## AIC: 177.08
##
## Number of Fisher Scoring iterations: 9
prob.glm <- predict(fit.glm,data_test,type = "response")
pred.glm <- rep(0,length(prob.glm))
pred.glm[prob.glm > 0.5] <- 1
table(pred.glm,data_test$crim01)
##
## pred.glm 0 1
## 0 65 9
## 1 9 69
mean(pred.glm != data_test$crim01)
## [1] 0.1184211
con los resultados obtenidos la tasa de error es de 13.15%
Con variables predictora con correlacion fuerte
fit.glm1 <- glm(crim01~indus+nox+age+dis+rad+tax+black+lstat,family = binomial,data_train)
summary(fit.glm1)
##
## Call:
## glm(formula = crim01 ~ indus + nox + age + dis + rad + tax +
## black + lstat, family = binomial, data = data_train)
##
## Deviance Residuals:
## Min 1Q Median 3Q Max
## -1.96371 -0.22572 -0.01948 0.00597 2.91220
##
## Coefficients:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -23.472808 5.213152 -4.503 6.71e-06 ***
## indus -0.052668 0.054397 -0.968 0.33293
## nox 45.719887 8.902831 5.135 2.81e-07 ***
## age 0.024557 0.012887 1.906 0.05671 .
## dis 0.265804 0.187558 1.417 0.15643
## rad 0.620643 0.146227 4.244 2.19e-05 ***
## tax -0.009174 0.003001 -3.057 0.00223 **
## black -0.007900 0.005400 -1.463 0.14348
## lstat -0.042885 0.040838 -1.050 0.29366
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for binomial family taken to be 1)
##
## Null deviance: 490.70 on 353 degrees of freedom
## Residual deviance: 160.24 on 345 degrees of freedom
## AIC: 178.24
##
## Number of Fisher Scoring iterations: 8
prob.glm1 <- predict(fit.glm1,data_test,type = "response")
pred.glm1 <- rep(0,length(prob.glm1))
pred.glm1[prob.glm1 > 0.5] <- 1
table(pred.glm1,data_test$crim01)
##
## pred.glm1 0 1
## 0 67 12
## 1 7 66
mean(pred.glm1 != data_test$crim01)
## [1] 0.125
con estos variables predictoras la tasa de error es de 15.13%.
LDA:
fit.lda <- lda(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat , data = data_train)
fit.lda
## Call:
## lda(crim01 ~ zn + indus + chas + nox + rm + age + dis + rad +
## tax + ptratio + black + lstat, data = data_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5056497 0.4943503
##
## Group means:
## zn indus chas nox rm age dis
## 0 22.460894 7.143631 0.04469274 0.4708587 6.399620 50.44469 5.153022
## 1 1.268571 15.137771 0.09142857 0.6337086 6.170234 85.79200 2.523521
## rad tax ptratio black lstat
## 0 4.201117 311.8156 17.92514 388.2874 9.549944
## 1 14.777143 506.8000 18.96800 319.7955 15.626343
##
## Coefficients of linear discriminants:
## LD1
## zn -0.0034199390
## indus 0.0102997449
## chas 0.0808703228
## nox 7.2261112533
## rm 0.0612382181
## age 0.0167656971
## dis -0.0027187266
## rad 0.0970467645
## tax -0.0018560309
## ptratio -0.0064709707
## black -0.0005065647
## lstat -0.0255231030
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$crim01)
##
## 0 1
## 0 69 15
## 1 5 63
mean(pred.lda$class != data_test$crim01)
## [1] 0.1315789
Con los resultados obtenidos con este metodo, la tasa de error es de 17.10% mas que la regresión logistíca.
Ahora con las variables predictoras con correlación fuerte.
fit.lda1 <- lda(crim01~indus+nox+age+dis+rad+tax+black+lstat , data = data_train)
fit.lda1
## Call:
## lda(crim01 ~ indus + nox + age + dis + rad + tax + black + lstat,
## data = data_train)
##
## Prior probabilities of groups:
## 0 1
## 0.5056497 0.4943503
##
## Group means:
## indus nox age dis rad tax black
## 0 7.143631 0.4708587 50.44469 5.153022 4.201117 311.8156 388.2874
## 1 15.137771 0.6337086 85.79200 2.523521 14.777143 506.8000 319.7955
## lstat
## 0 9.549944
## 1 15.626343
##
## Coefficients of linear discriminants:
## LD1
## indus 0.0118027470
## nox 7.2267222513
## age 0.0174936477
## dis -0.0237077126
## rad 0.0998803157
## tax -0.0021239136
## black -0.0005370087
## lstat -0.0291598636
pred.lda <- predict(fit.lda,data_test)
table(pred.lda$class,data_test$crim01)
##
## 0 1
## 0 69 15
## 1 5 63
mean(pred.lda$class != data_test$crim01)
## [1] 0.1315789
Con estas covariables se obtiene una tasa de error de 17.10% igual que la anterior pero por pricipio de parcimonia es mejor este modelo con pocas variables
KNN
set.seed(0511)
data_train$crim01 <- as.factor(data_train$crim01)
ctrl <- trainControl(method="repeatedcv",repeats = 3)
knnFit <- train(crim01~zn+indus+chas+nox+rm+age+dis+rad+tax+ptratio+black+lstat, data =data_train, method = "knn", trControl = ctrl, preProcess = c("center","scale"), tuneLength = 20)
knnFit
## k-Nearest Neighbors
##
## 354 samples
## 12 predictor
## 2 classes: '0', '1'
##
## Pre-processing: centered (12), scaled (12)
## Resampling: Cross-Validated (10 fold, repeated 3 times)
## Summary of sample sizes: 318, 319, 319, 319, 318, 319, ...
## Resampling results across tuning parameters:
##
## k Accuracy Kappa
## 5 0.9057672 0.8113559
## 7 0.8755556 0.7509298
## 9 0.8576984 0.7151367
## 11 0.8641799 0.7281249
## 13 0.8500000 0.6996443
## 15 0.8396296 0.6788571
## 17 0.8302646 0.6600574
## 19 0.8358201 0.6711574
## 21 0.8320635 0.6635468
## 23 0.8284127 0.6562289
## 25 0.8181746 0.6357683
## 27 0.8144444 0.6283471
## 29 0.8152116 0.6296212
## 31 0.8113492 0.6218266
## 33 0.8114286 0.6218624
## 35 0.8096032 0.6181100
## 37 0.8076984 0.6142139
## 39 0.7888889 0.5764674
## 41 0.7907672 0.5801909
## 43 0.7935714 0.5857363
##
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was k = 5.
Un K optimo para este metodo de clasificación es 5.
train.x <- as.matrix(data_train[,-c(1,14)])
test.x <- as.matrix(data_test[,-c(1,14)])
train.y <- as.matrix(data_train$crim01)
test.y <- as.matrix(data_test$crim01)
set.seed(0511)
pred.knn <- knn(train = train.x , test = test.x,cl = train.y,k=5)
table(pred.knn,test.y)
## test.y
## pred.knn 0 1
## 0 66 7
## 1 8 71
mean(pred.knn !=test.y)
## [1] 0.09868421
Este metodo tiene una tasa de error de 7.2% mucho mejor que los anteriores.
###Punto 7
library(MASS)
library(randomForest)
## randomForest 4.6-14
## Type rfNews() to see new features/changes/bug fixes.
##
## Attaching package: 'randomForest'
## The following object is masked from 'package:ggplot2':
##
## margin
train <- sample(1:nrow(Boston),length(Boston$crim)*0.70,replace = FALSE)
data.train <- Boston[train,-14]
data.test <- Boston[-train,-14]
y.train <- Boston[train,14]
y.test <- Boston[-train,14]
set.seed(0511)
boston.1 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
ntree = 500,mtry = ncol(Boston)-1)
boston.2 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
ntree = 500,mtry = (ncol(Boston)-1)/2)
boston.3 <- randomForest(x = data.train,y = y.train,xtest = data.test,ytest = y.test,
ntree = 500,mtry = sqrt(ncol(Boston)-1))
val <- data.frame(arboles = 1:500)
val["MSE1"] <- boston.1$test$mse
val["MSE2"] <- boston.2$test$mse
val["MSE3"] <- boston.3$test$mse
library(ggplot2)
library(reshape2)
xy <- melt(val,id=c("arboles"))
theme_set(theme_bw())
ggplot(xy)+geom_line(aes(x=arboles,y = value,color=variable))+
scale_color_manual(name=" m = ",labels=c("p","p/2","\u221Ap"),values = c("orange","red","green"))+labs(x = "Número de árboles",
y = "Error de clasificación") +
theme(legend.position="top")
Con esta grafica podemos observar que la prueba MSE es muy alta para un solo arbol, a mededida que aunmenta los arboles el error se estabiliza aproximademente cuando hay 150 arboles, La estaibilización se presenta para los tres valores de “m”. El m que menor reprensenta error es m = \(sqrt(p)\).Ojo que al variar la semilla se podra concluir diferente.
###Punto 8
(a) Divida el conjunto de datos en un conjunto de entrenamiento y un conjunto de prueba.
library(ISLR)
set.seed(0511)
train <- sample(1:nrow(Carseats),nrow(Carseats)*0.65,replace = FALSE)
carseats.train <- Carseats[train,]
carseats.test <- Carseats[-train,]
(b) Ajuste un árbol de regresión al conjunto de entrenamiento. Trace el árbol e interprete los resultados. ¿Qué prueba MSE obtienes?
library(rpart)
library(rpart.plot)
carseats.tree <- rpart(Sales~.,data = carseats.train)
rpart.plot(carseats.tree,box.palette = "RdBu",nn=TRUE)
pred <- predict(carseats.tree,carseats.test)
mean((pred - carseats.test$Sales)^2)
## [1] 5.463956
con el anterior resultado concluimos que el MSE es aproximadamente de 5.46.
(c) Utilice la validación cruzada para determinar el nivel óptimo de complejidad del árbol. ¿La poda del árbol mejora la prueba MSE?
printcp(carseats.tree)
##
## Regression tree:
## rpart(formula = Sales ~ ., data = carseats.train)
##
## Variables actually used in tree construction:
## [1] Advertising Age CompPrice Income Price ShelveLoc
## [7] US
##
## Root node error: 2084.7/260 = 8.0181
##
## n= 260
##
## CP nsplit rel error xerror xstd
## 1 0.328320 0 1.00000 1.00636 0.083418
## 2 0.118576 1 0.67168 0.68689 0.056595
## 3 0.045810 2 0.55310 0.61092 0.049569
## 4 0.041761 3 0.50729 0.59117 0.047620
## 5 0.030918 4 0.46553 0.56251 0.044079
## 6 0.030731 5 0.43461 0.55924 0.045740
## 7 0.023733 6 0.40388 0.53676 0.044968
## 8 0.021595 7 0.38015 0.53695 0.043152
## 9 0.019798 8 0.35856 0.53188 0.043446
## 10 0.018870 9 0.33876 0.53251 0.042926
## 11 0.015480 10 0.31989 0.51882 0.042220
## 12 0.015420 11 0.30441 0.51052 0.041593
## 13 0.013580 12 0.28899 0.50808 0.041005
## 14 0.012335 13 0.27541 0.50372 0.041579
## 15 0.010000 14 0.26307 0.50122 0.041944
plotcp(carseats.tree)
min <- which.min(carseats.tree$cptable[,"xerror"])
print(min)
## 15
## 15
cpt <- carseats.tree$cptable[which.min(carseats.tree$cptable[,"xerror"]),"CP"]
cpt
## [1] 0.01
con la validación cruzada obtenemos un tamaño de arbol igual a 13 y un cp de 0.01357961
ptree <- prune(carseats.tree,cp = cpt)
rpart.plot(ptree,box.palette = "RdBu",nn=TRUE)
pred1 <- predict(ptree,newdata = carseats.test)
m2 <- mean((pred1-carseats.test$Sales)^2)
m2
## [1] 5.463956
Se disminuye el MSE despues de que se poda.
(d) Utilice el enfoque de embolsado para analizar estos datos. ¿Qué prueba MSE obtienes? Use la función importance() para determinar qué variables son más importantes.
bag.carseats <- randomForest(Sales~.,carseats.train,mtry=(ncol(Carseats)-1) ,importance=TRUE)
pred2 <- predict(bag.carseats,newdata = carseats.test)
mean((pred2- carseats.test$Sales)^2)
## [1] 3.361658
notablemente mejora el MSE a 3.37627.
importance(bag.carseats)
## %IncMSE IncNodePurity
## CompPrice 26.67871283 183.54183
## Income 6.23358916 98.64888
## Advertising 13.57834499 110.41853
## Population 0.32338488 75.52763
## Price 68.47503499 599.01499
## ShelveLoc 77.26581305 725.54524
## Age 17.48598332 158.80999
## Education 0.01616168 57.99921
## Urban -1.71716623 11.09547
## US 5.07615837 11.71181
Se puede concluir con lo anterior que Price,Shelveloc,CompPrice son variables importantes para este analisis.
(e) Utilice bosques aleatorios para analizar estos datos. ¿Qué prueba MSE obtienes? Use la función importance() para determinar qué variables son las más importantes. Describa el efecto de m, el número de variables consideradas en cada división, sobre la tasa de error obtenida.
bag.carseats <- randomForest(Sales ~ ., data = carseats.train, mtry = sqrt((ncol(Carseats) - 1)), importance = TRUE)
pred.bag <- predict(bag.carseats, newdata = carseats.test)
mean((pred.bag - carseats.test$Sales)^2)
## [1] 3.391034
Con lo anterior tenemos un MSE de 3.38
importance(bag.carseats)
## %IncMSE IncNodePurity
## CompPrice 14.5403047 178.84227
## Income 2.6141067 154.02152
## Advertising 11.0107597 147.72102
## Population -0.9294335 133.60346
## Price 40.3920222 478.45992
## ShelveLoc 50.2416174 559.76511
## Age 10.8775741 194.24705
## Education -0.6566648 90.06653
## Urban -0.4877354 18.95984
## US 3.6856553 23.97304
Con lo anterior cuncluimos que las mejores varibales predictoras para el analisis son las mismas que el numeral anterior.
###Punto 9
(a) Cree un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones, y un conjunto de prueba que contenga las observaciones restantes.
train <- sample(x = 1:nrow(OJ),size = 800)
OJ.train <- OJ[train,]
OJ.test <- OJ[-train,]
(b) Ajuste un árbol a los datos de entrenamiento, con Compra como respuesta y las otras variables como predictores. Use la función summary() para generar estadísticas resumidas sobre el árbol y describa los resultados obtenidos. ¿Cuál es la tasa de error de entrenamiento? ¿Cuántos nodos terminales tiene el árbol?
set.seed(0511)
library(tree)
tree.oj <- tree(Purchase~.,data = OJ.train)
summary(tree.oj)
##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "SpecialCH" "WeekofPurchase"
## [5] "ListPriceDiff" "PctDiscMM"
## Number of terminal nodes: 10
## Residual mean deviance: 0.7071 = 558.6 / 790
## Misclassification error rate: 0.1525 = 122 / 800
Este arbalo tiene 8 nodos terminales y una tasa de error de entrenamiento de 0.1612
(c) Escriba el nombre del objeto del árbol para obtener una salida de texto detallada. Elija uno de los nodos terminales e interprete la información que se muestra.
tree.oj
## node), split, n, deviance, yval, (yprob)
## * denotes terminal node
##
## 1) root 800 1063.000 CH ( 0.61875 0.38125 )
## 2) LoyalCH < 0.5036 351 424.900 MM ( 0.29345 0.70655 )
## 4) LoyalCH < 0.0356415 47 0.000 MM ( 0.00000 1.00000 ) *
## 5) LoyalCH > 0.0356415 304 389.300 MM ( 0.33882 0.66118 )
## 10) SalePriceMM < 2.04 167 176.600 MM ( 0.22156 0.77844 )
## 20) SpecialCH < 0.5 137 124.000 MM ( 0.16788 0.83212 ) *
## 21) SpecialCH > 0.5 30 41.460 MM ( 0.46667 0.53333 ) *
## 11) SalePriceMM > 2.04 137 189.700 MM ( 0.48175 0.51825 )
## 22) LoyalCH < 0.168934 27 22.650 MM ( 0.14815 0.85185 ) *
## 23) LoyalCH > 0.168934 110 150.700 CH ( 0.56364 0.43636 )
## 46) WeekofPurchase < 244.5 17 15.840 MM ( 0.17647 0.82353 ) *
## 47) WeekofPurchase > 244.5 93 122.100 CH ( 0.63441 0.36559 ) *
## 3) LoyalCH > 0.5036 449 341.700 CH ( 0.87305 0.12695 )
## 6) LoyalCH < 0.764572 193 220.800 CH ( 0.74093 0.25907 )
## 12) ListPriceDiff < 0.235 74 102.600 MM ( 0.50000 0.50000 )
## 24) PctDiscMM < 0.196196 59 78.900 CH ( 0.61017 0.38983 ) *
## 25) PctDiscMM > 0.196196 15 7.348 MM ( 0.06667 0.93333 ) *
## 13) ListPriceDiff > 0.235 119 82.090 CH ( 0.89076 0.10924 ) *
## 7) LoyalCH > 0.764572 256 64.200 CH ( 0.97266 0.02734 ) *
Escogemos el nodo terminal 15, debido a que es un nodo terminal debido al asterisco.apreciamos que su criterio de división es LoyalCH > 0.764572 esta tiene 254 observaciones en este nodo y una desviación estandar 64.090. La predicción general para la rama CH,también podemos ver que el 97.24% de las observaciones toman el valor de CH, mientras que el 0.02756 tomal el valor de MM.
(d) Cree un diagrama del árbol e interprete los resultados.
plot(tree.oj)
text(tree.oj, pretty = 0)
Podemos apreciar que el nodo principal es LoyalCH si es menor o igual a 0.50 solo usa a su misma variables , pero si es menor utiliza a PriceDiff.
(e) Predecir la respuesta en los datos de prueba, y producir una matriz de confusión comparando las etiquetas de prueba con las etiquetas de prueba predichas. ¿Cuál es la tasa de error de prueba?
library(caret)
tree.pred <- predict(tree.oj, OJ.test, type = "class")
matrix<-confusionMatrix(tree.pred, OJ.test$Purchase)
matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 138 33
## MM 20 79
##
## Accuracy : 0.8037
## 95% CI : (0.7512, 0.8494)
## No Information Rate : 0.5852
## P-Value [Acc > NIR] : 1.923e-14
##
## Kappa : 0.5887
##
## Mcnemar's Test P-Value : 0.09929
##
## Sensitivity : 0.8734
## Specificity : 0.7054
## Pos Pred Value : 0.8070
## Neg Pred Value : 0.7980
## Prevalence : 0.5852
## Detection Rate : 0.5111
## Detection Prevalence : 0.6333
## Balanced Accuracy : 0.7894
##
## 'Positive' Class : CH
##
Observamos que la precisión del modelo es de 82.89%, siendo mas preciso para predecir a CH con un 90.65% mientras que para predecir MM sólo es de un 70.91%. (f) Aplique la función cv.tree() al conjunto de entrenamiento para determinar el tamaño óptimo del árbol.
cv.oj <- cv.tree(tree.oj,FUN = prune.misclass)
cv.oj
## $size
## [1] 10 9 6 2 1
##
## $dev
## [1] 142 144 151 157 305
##
## $k
## [1] -Inf 0.000000 4.333333 6.250000 145.000000
##
## $method
## [1] "misclass"
##
## attr(,"class")
## [1] "prune" "tree.sequence"
(g) Produzca un diagrama con el tamaño del árbol en el eje xy la tasa de error de clasificación con validación cruzada en el eje y.
library(ggplot2)
library(reshape2)
valores <- data.frame(Size=cv.oj$size,Dev=cv.oj$dev)
theme_set(theme_bw())
ggplot(data = valores, aes(x=Size, y=Dev)) +
geom_line(color="red", linetype = "dashed") +
geom_point() +
scale_x_discrete(limits=c(1:9))
(h) ¿Qué tamaño de árbol corresponde a la tasa de error de clasificación con validación cruzada más baja? 5 es el tamaño con menor error
(i) Produzca un árbol podado correspondiente al tamaño óptimo del árbol obtenido mediante validación cruzada. Si la validación cruzada no conduce a la selección de un árbol podado, cree un árbol podado con cinco nodos terminales.
prune.oj <- prune.misclass(tree.oj,best = 5)
plot(prune.oj)
text(prune.oj,pretty = 0)
(j) Compare las tasas de error de entrenamiento entre los árboles podados y no podados. ¿Cuál es más alto?
summary(prune.oj)
##
## Classification tree:
## snip.tree(tree = tree.oj, nodes = c(10L, 3L))
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "WeekofPurchase"
## Number of terminal nodes: 6
## Residual mean deviance: 0.8552 = 679 / 794
## Misclassification error rate: 0.1688 = 135 / 800
summary(tree.oj)
##
## Classification tree:
## tree(formula = Purchase ~ ., data = OJ.train)
## Variables actually used in tree construction:
## [1] "LoyalCH" "SalePriceMM" "SpecialCH" "WeekofPurchase"
## [5] "ListPriceDiff" "PctDiscMM"
## Number of terminal nodes: 10
## Residual mean deviance: 0.7071 = 558.6 / 790
## Misclassification error rate: 0.1525 = 122 / 800
observamos que hay una diferencia en el erro mayor el arbol podado. (k) Compare las tasas de error de prueba entre los árboles podados y no podados. ¿Cuál es más alto?
prune.pred <- predict(prune.oj, OJ.test, type = "class")
matrix<-confusionMatrix(prune.pred, OJ.test$Purchase)
matrix
## Confusion Matrix and Statistics
##
## Reference
## Prediction CH MM
## CH 140 39
## MM 18 73
##
## Accuracy : 0.7889
## 95% CI : (0.7353, 0.836)
## No Information Rate : 0.5852
## P-Value [Acc > NIR] : 1.161e-12
##
## Kappa : 0.553
##
## Mcnemar's Test P-Value : 0.008071
##
## Sensitivity : 0.8861
## Specificity : 0.6518
## Pos Pred Value : 0.7821
## Neg Pred Value : 0.8022
## Prevalence : 0.5852
## Detection Rate : 0.5185
## Detection Prevalence : 0.6630
## Balanced Accuracy : 0.7689
##
## 'Positive' Class : CH
##
observamos que bajo la presicion del modelo notablemente con respecto al árbol sin podar.
###Punto 10
(a) Elimine las observaciones para las que se desconoce la información salarial, y luego transforme logarítmicamente los salarios.
Hitters <- na.omit(Hitters)
Hitters$Salary <- log(Hitters$Salary)
(b) Cree un conjunto de entrenamiento que consista en las primeras 200 observaciones, y un conjunto de prueba que consista en las observaciones restantes.
train <- 1:200
Hitters.train <- Hitters[train,]
Hitters.test <- Hitters[-train,]
(c) Realice un refuerzo en el conjunto de entrenamiento con 1,000 árboles para un rango de valores del parámetro de contracción \(\lambda\). Produzca un gráfico con diferentes valores de contracción en el eje xy el conjunto de entrenamiento MSE correspondiente en el eje y.
library(gbm)
## Loaded gbm 2.1.5
set.seed(0511)
lambdas = seq(0.001, 0.3, 0.005)
mse <- rep(0, length(lambdas))
for (i in 1:length(lambdas)) {
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
pred.train <- predict(boost.hitters, Hitters.train, n.trees = 1000)
mse[i] <- mean((pred.train - Hitters.train$Salary)^2)
}
library(ggplot2)
mse_graph <- data.frame(lambdas)
mse_graph["MSE"]<- mse
names(mse_graph) <- c("Lambdas", "MSE")
theme_set(theme_bw())
ggplot(data = mse_graph, aes(x = Lambdas, y = MSE)) +
geom_line(color="red", linetype = "dashed") +
geom_point()
(d) Produzca un gráfico con diferentes valores de contracción en el eje xy el conjunto de prueba MSE correspondiente en el eje y.
mse2 <- rep(NA, length(lambdas))
for (i in 1:length(lambdas)) {
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[i])
yhat <- predict(boost.hitters, Hitters.test, n.trees = 1000)
mse2[i] <- mean((yhat - Hitters.test$Salary)^2)
}
library(ggplot2)
mse_graph2 <- data.frame(lambdas)
mse_graph2["MSE"]<- mse2
names(mse_graph2) <- c("Lambdas", "MSE")
theme_set(theme_bw())
ggplot(data = mse_graph2, aes(x = Lambdas, y = MSE)) +
geom_line(color="red", linetype = "dashed") +
geom_point()
min(mse2)
## [1] 0.2461295
lambdas[which.min(mse2)]
## [1] 0.131
(e) Compare la prueba MSE de refuerzo con la prueba MSE que resulta de aplicar dos de los enfoques de regresión vistos en los Capítulos 3 y 6.
mod8.10 <- lm(Salary~.,data = Hitters.train)
pred8.10 <- predict(mod8.10,Hitters.test)
mean((pred8.10-Hitters.test$Salary)^2)
## [1] 0.4917959
library(pls)
##
## Attaching package: 'pls'
## The following object is masked from 'package:caret':
##
## R2
## The following object is masked from 'package:stats':
##
## loadings
pcr.fit=pcr(Salary~., data=Hitters.train , scale=TRUE , validation ="CV")
validationplot(pcr.fit ,val.type="MSEP")
pcr.pred=predict (pcr.fit ,Hitters.test,ncomp =1)
mean((pcr.pred -Hitters.test$Salary)^2)
## [1] 0.4661183
observamos que el MSE de la regresion lineal y de las componentes principales es mayor el MSE por el metodo Boosting
(f) ¿Qué variables parecen ser los predictores más importantes en el modelo impulsado?
boost.hitters <- gbm(Salary ~ ., data = Hitters.train, distribution = "gaussian", n.trees = 1000, shrinkage = lambdas[which.min(mse2)])
summary(boost.hitters)
## var rel.inf
## CAtBat CAtBat 22.1891184
## PutOuts PutOuts 8.1156374
## CRBI CRBI 7.7992242
## CWalks CWalks 7.2267797
## Walks Walks 6.9703139
## Assists Assists 6.0854321
## CHmRun CHmRun 5.2576243
## Hits Hits 5.2422222
## Years Years 4.8380328
## CRuns CRuns 4.6471959
## RBI RBI 4.6430518
## AtBat AtBat 3.5618987
## HmRun HmRun 3.5243625
## Errors Errors 3.0632824
## CHits CHits 2.7467221
## Runs Runs 2.7155911
## Division Division 0.6010006
## NewLeague NewLeague 0.5563958
## League League 0.2161140
En la tabla se observa que la variable mas importante es CatBat.
(g) Ahora aplique el embolsado al conjunto de entrenamiento. ¿Cuál es el conjunto de prueba MSE para este enfoque?
bag.hitters <- randomForest(Salary ~ ., data = Hitters.train, mtry = 19)
yhat.bag <- predict(bag.hitters, newdata = Hitters.test)
mean((yhat.bag - Hitters.test$Salary)^2)
## [1] 0.228878
El MSE es de 0.2339 es mucho mejor que los otros metodos por el momento. ###Punto 11
(a) Cree un conjunto de entrenamiento que consista en las primeras 1,000 observaciones,y un conjunto de prueba que consta de las observaciones restantes.
(b) Ajuste un modelo de refuerzo al conjunto de entrenamiento con Compra como respuesta y las otras variables como predictores. Use 1,000 árboles y un valor de contracción de 0.01. ¿Qué predictores parecen ser los más importantes?
(c) Utilice el modelo de refuerzo para predecir la respuesta en los datos de la prueba. Predecir que una persona realizará una compra si la probabilidad estimada de compra es superior al 20%. Forme una matriz de confusión. ¿Qué fracción de las personas que predijeron hacer una compra en realidad la hacen? ¿Cómo se compara esto con los resultados obtenidos al aplicar KNN o regresión logística a este conjunto de datos?
###Punto 12
###Punto 4
library(e1071)
set.seed(1)
x <- rnorm(100)
y <- 4 * x^3 + 1 + rnorm(100)
class <- sample(100, 50)
y[class] <- y[class] + 3
y[-class] <- y[-class] - 3
plot(x[class], y[class], col = "green", xlab = "X", ylab = "Y", ylim = c(-10, 25), xlim = c(-1.5, 2.0))
points(x[-class], y[-class], col = "blue")
z <- rep(-1, 100)
z[class] <- 1
data <- data.frame(x = x, y = y, z = as.factor(z))
train <- sample(100, 50)
data.train <- data[train, ]
data.test <- data[-train, ]
svm.linear <- svm(z ~ ., data = data.train, kernel = "linear", cost = 10)
plot(svm.linear, data.train)
table(predict = predict(svm.linear, data.train), truth = data.train$z)
## truth
## predict -1 1
## -1 21 1
## 1 4 24
svm.poly <- svm(z ~ ., data = data.train, kernel = "polynomial", cost = 10)
plot(svm.poly, data.train)
table(predict = predict(svm.poly, data.train), truth = data.train$z)
## truth
## predict -1 1
## -1 23 1
## 1 2 24
svm.radial <- svm(z ~ ., data = data.train, kernel = "radial", gamma = 1, cost = 10)
plot(svm.radial, data.train)
table(predict = predict(svm.radial, data.train), truth = data.train$z)
## truth
## predict -1 1
## -1 25 0
## 1 0 25
plot(svm.linear, data.test)
table(predict = predict(svm.linear, data.test), truth = data.test$z)
## truth
## predict -1 1
## -1 24 0
## 1 1 25
plot(svm.poly, data.test)
table(predict = predict(svm.poly, data.test), truth = data.test$z)
## truth
## predict -1 1
## -1 23 1
## 1 2 24
plot(svm.radial, data.test)
table(predict = predict(svm.radial, data.test), truth = data.test$z)
## truth
## predict -1 1
## -1 19 0
## 1 6 25
###Punto 5
(a) Genere un conjunto de datos con n = 500 y p = 2, de modo que las observaciones pertenezcan a dos clases con un límite de decisión cuadrático entre ellas. Por ejemplo, puede hacer esto de la siguiente manera:
> x1 = runif (500) -0.5
> x2 = runif (500) -0.5
> y = 1 * (x1 ^ 2-x2 ^ 2> 0)
(b) Grafique las observaciones, coloreadas de acuerdo con sus etiquetas de clase. Su diagrama debe mostrar X1 en el eje xy X2 en el eje y.
(c) Ajuste un modelo de regresión logística a los datos, utilizando X1 y X2 como predictores.
(d) Aplique este modelo a los datos de entrenamiento para obtener una etiqueta de clase predicha para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas. El límite de decisión debe ser lineal.
(e) Ahora ajuste un modelo de regresión logística a los datos utilizando funciones no lineales de X1 y X2 como predictores (por ejemplo, X12, X1 × X2, log (X2), etc.).
(f) Aplique este modelo a los datos de entrenamiento para obtener una etiqueta de clase predicha para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas. El límite de decisión debe ser obviamente no lineal. Si no es así, repita (a) - (e) hasta que encuentre un ejemplo en el que las etiquetas de clase predichas sean obviamente no lineales.
(g) Ajuste un clasificador de vector de soporte a los datos con X1 y X2 como predictores. Obtenga una predicción de clase para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas.
(h) Ajuste un SVM usando un núcleo no lineal a los datos. Obtenga una predicción de clase para cada observación de entrenamiento. Trace las observaciones, coloreadas de acuerdo con las etiquetas de clase predichas.
(i) Comente sus resultados.
###Punto 6
(a) Genere datos de dos clases con p = 2 de tal manera que las clases sean apenas separables linealmente.
(b) Calcule las tasas de error de validación cruzada para los clasificadores de vectores de soporte con un rango de valores de costo. ¿Cuántos errores de capacitación están mal clasificados para cada valor de costo considerado, y cómo se relaciona esto con los errores de validación cruzada obtenidos?
(c) Genere un conjunto de datos de prueba apropiado y calcule los errores de prueba correspondientes a cada uno de los valores de costo considerados. ¿Qué valor del costo conduce a la menor cantidad de errores de prueba, y cómo se compara esto con los valores de costo que producen la menor cantidad de errores de capacitación y la menor cantidad de errores de validación cruzada?
(d) Discuta sus resultados.
###Punto 7
(a) Cree una variable binaria que tome un 1 para los automóviles con millaje de gasolina por encima de la mediana, y un 0 para automóviles con millaje de gasolina por debajo de la mediana.
(b) Ajuste un clasificador de vector de soporte a los datos con varios valores de costo, para predecir si un automóvil obtiene un millaje de gasolina alto o bajo. Informe los errores de validación cruzada asociados con diferentes valores de este parámetro. Comenta tus resultados.
(c) Ahora repita (b), esta vez usando SVM con núcleos de base radial y polinomial, con diferentes valores de gamma y grado y costo. Comenta tus resultados.
(d) Haga algunos gráficos para respaldar sus afirmaciones en (b) y (c). Sugerencia: En el laboratorio, utilizamos la función plot() para objetos svm solo en casos con p = 2. Cuando p> 2, puede usar la función plot() para crear gráficos que muestren pares de variables a la vez. Esencialmente, en lugar de escribir
> plot (svmfit, dat)
**donde svmfit contiene su modelo ajustado y dat es un marco de datos que contiene sus datos, puede escribir*
> plot (svmfit, dat, x1∼x4)
para graficar solo las variables primera y cuarta. Sin embargo, debe reemplazar x1 y x4 con los nombres correctos de las variables. Para obtener más información, escriba ?Plot.svm.
###Punto 8
(a) Cree un conjunto de entrenamiento que contenga una muestra aleatoria de 800 observaciones, y un conjunto de prueba que contenga las observaciones restantes.
(b) Ajuste un clasificador de vector de soporte a los datos de entrenamiento usando cost = 0.01, con Compra como respuesta y las otras variables como predictores. Use la función summary() para generar estadísticas resumidas y describir los resultados obtenidos.
(c) ¿Cuáles son las tasas de error de capacitación y prueba?
(d) Use la función tune() para seleccionar un costo óptimo. Considere values en el rango de 0.01 a 10.
(e) Calcule las tasas de error de entrenamiento y prueba usando este nuevo valor por el cost
(f) Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo radial Use el valor predeterminado para gamma.
(g) Repita las partes (b) a (e) usando una máquina de vectores de soporte con un núcleo polinomial Establecer grado = 2.
(h) En general, ¿qué enfoque parece dar los mejores resultados con estos datos?